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Abstract 

This paper introduces a novel variational approach for image compression motivated 
by recent PDE-based approaches combining edge detection and Laplacian inpainting. The 
essential feature is to encode the image via a sparse vector field, ideally concentrating on 
a set of measure zero. An equivalent reformulation of the compression approach leads to a 
variational model resembling the ROF-model for image denoising, hence we further study 
the properties of the effective regularization functional introduced by the novel approach 
and discuss similarities to TV and TGV functionals. Moreover we computationally inves¬ 
tigate the behaviour of the model with sparse vector fields for compression in particular 
for high resolution images and give an outlook towards denoising. 

Keywords: Image compression, denoising, reconstruction, diffusion inpainting, sparsity, 
total variation 


1 Introduction 

Image compression is a topic of interest since the beginning of digital imaging, remaining 
relevant continuously due to the ongoing improvement in image resolution. While standard 
approaches are based on orthogonal bases and frames like cosine transforms or wavelets, an 
alternative route based on ideas from partial differential equations has emerged recently (cf. 
naan). In the latter case particular attention is paid to compressions from which cartoons 
can be reconstructed accurately avoiding the artefacts of the above mentioned standard ap¬ 
proaches by a direct treatment of edges. Roughly speaking, their idea is to detect edges and 
store the image value in pixels on both sides of an edge. The remaining parts of the image 
are completed by harmonic inpainting. An alternative interpretation of the two-sided pixel 
values, worked out more clearly in the osmosis setting of [20] . is that a vector field v corre¬ 
sponding to the normal derivatives of the image at the edge location is stored. The inpainting 
of the image u then corresponds to a solution of 

A u = V • v in O , (1) 

where ft is the image domain into which v is extended by zero off the edges. 
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In this paper we start from reinterpretation of the PDE-based compression in terms of 
sparsity, which directly translates into a variational framework. The key observation is that 
with the edge detection and zero extension of v one essentially looks for a sparse vector field 
that leads to a certain precision in the reconstruction via 0. In the spirit of the predomi¬ 
nant sparsity regularization of imaging problems we study a direct variational approach: we 
minimize an lA-type norm of the vector field subject to the constraint that u reconstructed 
via 0 approximates a given image / up to a certain tolerance, i.e. 


ll'L’Hi —» min subject to A u = V • v , \\u — f\\ 2 < e. (2) 

11, V 

As we shall discuss below a more rigorous statement of the problem takes into account that 
v needs to be interpreted as a vectorial Radon measure on Q (similar to gradients of BV- 
functions) in a continuum setting. The properties of a limiting continuum model appear 
to be of particular advantage for the compression issues, since one expects to concentrate 
the measure v on a set of Lebesgue measure zero. This means that for a suitable discrete 
approximation of the model the number of pixels we need to store the vector field in divided 
by the total number of pixels tends to zero. A direct consequence is an increase in the 
compression rates as the image resolution increases, a highly desirable property. 

Apart from compression we shall investigate models like the above one as a regularization 
for more general imaging problems. The relation to denoising models can readily be seen 
when the above constraint of approximating / is incorporated via a Lagrange functional. 
Indeed, there exists a Lagrange parameter A > 0 such that © is equivalent to 

~\\u — /|| 2 + ||v||i —y min subject to Au = V • v . (3) 

2 u,v 

Hence, we may also interpret the norm of v as an implicit regularization of v and simply 
replace the first term by an arbitrary data term to treat other imaging problems. This is 
more apparent when we replace the constraint by a natural special solution v = \7u. In this 
case ([3]) is just the ROF-model for denoising (cf. [T8] ) and indeed an equivalence relation 
holds in spatial dimension one. In higher spatial dimensions it becomes apparent that a key 
role is played by the curl V x v. If the curl of v vanishes, it simply becomes a gradient 
vector field and hence we again recover the ROF-model. This motivates to study further 
generalizations of the model also penalizing Vxu. With the interpretation that v is related 
to the gradient of u, the additional term becomes a higher order regularization effectively. 
Functionals of this type are currently studied in particular to reduce staircasing artefacts in 
total variation regularization (cf. [9, BJ E]), and indeed we shall be able to draw very close 
analogies to the recently very popular TGV-approach (cf. [3]). The reduction of staircasing 
is confirmed by computational experiments. However, in the denoising case we shall see that 
the divergence part alone does not suffice for appropriate smoothing. 

The remainder of the paper is organized as follows: In Section 2 we discuss the model and 
its variants. Then, we proceed to a discussion of some theoretical properties in Section 3, which 
provide some insights into the sparse vector field model. Section 4 introduces a numerical 
solution based on primal-dual optimization methods, which is used for some experimental 
studies in Section 5. Finally, we provide a conclusion and directions for future research. 
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2 Variational Model and Regularization Functionals 


In order to obtain an appropriate formulation of © we proceed as in [5j and interpret v as a 
d-dimensional Radon measure on 0 C M. d . The regularization functional is then 


’II.M(Q) = sup <p- dv, 

<peC(n) d ,\\ip\\oo<i Jq 


(4) 


where © is to be understood in a weak form as well. Hence, © is rewritten 


as 


^\\u — fllo + IMI aa(o) min 

2 " J " 2 11 ll>,(S2) ueL*(n),veM(ny 


subject to A u = V-v. 


(5) 


The model can be formulated as in recent approaches for denoising by defining w — \7u — v 
(in the sense of distributions). With xo being the characteristic function of the set {0} we 
obtain 


^\\u-f\\l + \\^u-w\\ M{n) + xo{^-w) ^ min 

2 v ' ueL 2 (n),weT>'(n) d 


( 6 ) 


One observes that the regularization functional is now an infimal convolution of total variation 
and a functional of Vw. The same structure is apparent in the recently popularized TGV- 
model (cf. [3]), which in the analogous setting reads 


A, 


1“ - /III + IIVu - HI M (Q) + l|VHU ( n) , “ in . 

uGL j (Q) ,weA4(Q) d 


(7) 


A major difference of the TGV approach to our new model is the fact that in our case only 
the divergence of v respectively w is penalized, which might be too weak to achieve suitable 
regularization properties. This can obviously be realized if we add additional regularization 
terms depending on V x v, which is natural since a divergence-free vector field is constant 
if and only if its curl vanishes. Note that V x \7u = 0, hence V x v = V x w, i.e. we can 
formulate regularization either on v or on w. The most general formulation is given by 

t;\\u- f\\l + \\Vu-w\\ Mm + F(V-w) + G(V xw) -> min (8) 

2 v ueL 2 (Q) i weT> f (Q) d 

with convex functionals F and G, which however exceeds the scope of this paper and is left 
as subject of future research. 

We finally recast the above results in terms of the regularization they induce on u. In 
this respect we also discuss the boundary conditions in 0. respectively its weak formulation. 
Natural boundary conditions are no-flux conditions (\7u — v) • n = 0 on 90, which means the 
used weak formulation of ([T]) is 


/ A cp udx+ V tp • dv(x) = 0 \/(p G C 2 (0), Vcp • n — 0 on 90 . 

Jq Jq 

Hence, we can define the regularization functional R : L 1 (O) [0, oo] 


(9) 


R(u) := inf ||w||jw(n) 

v satisfying (I9J 


(10) 
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Once we have defined the regularization terms it is straight-forward to extend the vari¬ 
ational model to other imaging tasks, e.g. by just changing the data fidelity. Moreover, we 
can consider Bregman iterations (cf. [16]) 

u-f\\l +R(u) - {p k ,u)\ , p k eR{u k ), (ll) 

as well as other scale space methods such as the gradient flow (cf. [1]) dtu G —dR(u) and the 
inverse scale space method (cf. |6]). 


U _Li . I A 

u G argmin [ — 


3 Properties of Regularization by Sparse Vector Fields 


In the following we further discuss some properties of the regularization functional R defined 
via (10). To avoid obvious technicalities with constants, we restrict ourselves to the space 


Ll(n) = { u e L\n) 


/ u dx = 0 } 

Jn 


( 12 ) 


if is a bounded domain. 

We start with some topological properties induced by R: 

Theorem 1. Let Q be a sufficiently regular domain. Then there exists a constant c > 0 such 
that 

IMIl^q) < cR(u ) (13) 

for all u € Ll(Q). Moreover, dom(R) is a subspace o/L^O) and R is a norm on dom(R) fl 
L\{ P). Finally, 

R(u) < TV(u ) Vue BV(Q). (14) 


Proof. We have 


Mki(n) = SU P 

4>eL°°(ci),\\4 


=<i 


/ u{x)cj){x) dx. 

Jn 


(15) 


For <f> G L°°(fl) we define w as the weak solution of the Poisson equation —Aw = <f> with 
homogeneous Neumann boundary conditions and mean value zero. Thus, using the weak 
formulations we have 

/ u(x)<f>(x) dx= \7w(x) • dv(x). (16) 

Jq Jq 

Regularity of solutions of the Poisson equation yields continuity of w and the existence of a 
constant c such that HVwHoo < cH^Hoo = c, for all <f> G L°°(Q). Hence, 




sup / u(x)(j>(x ) dx < cllull^m , 

fi),|Woo<Wfi 


(17) 


which yields the estimate of the L^norm. 

The one-homogeneity and triangle inequality follow in a straight-forward way from the 
definition, consequently R is a norm on a subspace of Estimate (14) is obtained since 

v = X7u satisfies Q, hence the infimum over all admissible v is less or equal to the total 
variation. □ 
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A next step towards the understanding of properties of R is an investigation of its sub¬ 
differential, with subsequent consequences for optimality conditions of For brevity we 
use a formal approach based on Lagrange multipliers. We have p G dR{u) if and only if 
p — d u L{u , v, q), for solutions (v, q ) of the saddle-point problem 


inf sup L(u, v , q) 
v q 


(18) 


for given u and the Lagrangian is defined as 


L{u,v,q) 



A q u dx + 


/ \7q • dv(x) . 

Jn 


(19) 


Thus, we find p = A q and the optimality conditions for the saddle point problem yield © 
and — Vg G dlMl^n). 

It is instructive to compare the subgradients of R with those of TV. Indeed with similar 
reasoning one can show that p G dTV(u) if p = —V • g for a vector field g G $|M|.m(0) and 
v = Vu. This means that if we can write g — —Vq we also obtain p G dR(u). In particular 
this opens the door towards a simple verification whether solutions of the ROF-model are 
also solutions of the sparse vector field model One simply has to inspect the subgradient 
in the optimality condition and check whether the associated vector field g can be written 
as a gradient. We will exemplify this in the case of the most well-known example for the 
ROF model, the reconstruction of the indicator function of a ball on D = M. d (cf. [EJ). This 
function is an eigenfunction of TV, i.e. there exists A (depending on the radius R of the ball) 
such that 

A U = V -ge dTV(u ). (20) 


It is easy to see that g 
and F satisfies 


VF(6), where b is the signed distance function of the ball (cf. [TO] ) 


F'(z) 



if z <0 
if z > 0. 


( 21 ) 


Hence, we have g = — X7q for q = — F(6), which implies that Xu G dR{u ), with the same 
value. The results in [3j imply that the variational model © reconstructs data / being the 
indicator function of a ball in the form u — cf , with c < 1 depending on A and a. Moreover, 
the Bregman iteration and inverse scale space methods reconstruct / exactly after a finite 
number of iterations respectively finite time. 

Finally we return to the original idea of compressing an image by encoding a sparse 
vector field. For this sake it is desireable that v has support on a set of small (or even 
zero) Lebesgue measure. Thinking about the continuum case as a limit of discrete pixel 
images, the asymptotic property of zero Lebesgue measure means that the image (in 2D) 
can be encoded by a number of values proportional to the square root of the number of 
pixels. Consequently the compression rate of such a PDE-based approach should improve 
with higher image resolution, which is highly relevant given the current trend of screen and 
camera resolution. We already see from the example of the indicator function of a ball 
above that we can expect the method to encode a piecewise constant image by vector fields 
concentrated on the edge sets. For more complicated images the vector field potentially needs 
to have a larger support to obtain a suitable reconstruction, since away from the support 
of v the function u is just harmonic. A better understanding of the compression properties 
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would need a characterization of the structure of minimizers, similar to [7, QjJj. While the 
one-dimensional case is equivalent to total variation regularization and hence always yields v 
concentrated on a set of zero Lebesgue measure (cf. [17]), the multi-dimensional case is less 
clear and left as an interesting topic of future research. Further studies on the compression 
properties will be carried out below by computational experiments. 


4 Numerical Solution via Primal-Dual Methods 


In order to solve our minimization problem ([5]) numerically, we at first need to discretize it. 
Thereto we will adopt the notation of the continuous functions it, v, f and the operators V, 
V- and A, but from now on, we are thereby referring to their discretized versions, which we 
shall comment on in the following. 

For simplification, we assume the normalized images to be quadratic, i.e. /, u G [0, l] ArxiV . 
The pixel grid can be written as {(i/q jh): 1 < i, j < V}, where h denotes the spacing size. 
We use forward finite differences with Neumann boundary conditions for the discretization of 
the gradient of u and in order to preserve the adjoint structure the divergence is discretized 
with backward finite differences. Moreover, in this discrete setting the d-dimensional Radon 
measure on ft C M. d becomes the discrete Ll-norm. Considering v as being related to the 
gradient of iq the regularization term in problem ^ can be interpreted as the discrete total 
variation norm. We decided to use its isotropic version which is given by 

iivuiii = Ekv^-i = E A v Ai) 2 + (( Vw )l) 2 - 

h3 ij 

Consequently, the minimization problem ^ reads: 

— \\u ~ f II 2 + IMIi min subject to A u = V • v 

2 U,v 

or equivalently: 

h|«-/||| + |M|i + Xo(Au -V-v) —> min, (24) 

Z u,v 

where xo is again the characteristic function of the set {0}. Defining 


( 22 ) 


(23) 


x := (iq v) rj 


G(x) := ^\\u - /||1 + 


|i, F(Kx ) := Xo(A u -V-v), 


(25) 


one can easily see that we can calculate a solution of the above problem (24) by applying a 
version of the recently very popular primal-dual algorithms (cf. [TTj l8]) designed for efficiently 
solving general minimization problems of the form 


F[Kx ) + G{x) min, (26) 

X 

where F and G are proper convex lower-semicontinuous functionals. 

We decided to use the first-order primal-dual algorithm as proposed by Chambolle and 
Pock (cf. [8]) given by Algorithm [lj Adopting their notation we can now derive the updates 
concerning our minimization problem ( [24] ). Thereto we at first have to calculate the dual 
functional F*(y ) = sup X {(y,x) — F(x)}. Since in our case F is the characteristic function 
of the set {0}, it is straightforward to see that F* equals zero. Hence the update for the 
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Algorithm 1 Primal-Dual Algorithm by Chambolle and Pock 
Input: r, a > 0, 6 G [0,1] 

Initialization: x°, y°, x° = x° 
for n > 0 do 

y n+1 = (J + adF*) -1 (y n + aKx n ) 
x n+1 = (J + r9G) _1 (x n - rFT*y n+1 ) 

^n+l = a.n+1 + Q ( x n +1 _ 

end for 


dual variable ?/ simplifies to 7/ n+1 — y n -\- aKx n . Next we consider the update for the primal 
variable x. As the subdifferentials of G with respect to u and v are independent of v and 
u , respectively, we can update each component of x separately (see Algorithm 2 ). Using the 
norm of the sum of the primal and dual residual given by (cf. |l 2 j) 

||^pn-|-l ||2 _|_ ii^n-l-11|2 (27) 


where 


pn +1 _ 


x n _ x n +1 


n _ n +1 

- iT(y n - y n+1 ), £> n+1 = ^- y - - K(x n - x n+1 ), (28) 

a 


as a stopping criterion the implementation of our minimization problem (24) can be sum¬ 
marized by Algorithm [ 2 | where the two-dimensional isotropic shrinkage operator is defined 
by: 


shrinkage(z, 7 ) = max(||;z||i — 7 )- 


\z\\i 


(29) 


Algorithm 2 Primal-Dual Algorithm for Minimization of (24) 


Input: image /, A > 0, r, a > 0, 6 G [0,1], max no of iterations, e > 0 
Initialization: v °, ?/ 0 , 4Z° = ?x°, v° = ?; 0 

while primal-dual residual > e and n < max no of iterations do 

y n+1 = (I + <7<9F*) _1 (y n + crV • (W 1 - v n )) = y n + (TV • (V« n - F 1 ) 
■u ra+1 = (/ + r9 u G) _1 (« n - rAy n+1 ) = ^ (Ar/ + u n - rAy n+1 ) 
v n+1 = (I + rd d G)~ l ( v n - rVy n+1 ) = shrinkage^" - rVy n+1 , r) 
y n+ i = u n+1 + 9 ( u n+1 - u n ) 

^n+l = v n +1 + y ^n+1 _ 

end while 


As already mentioned in Section [ 2 ] the model discussed so far can be further extended 
by considering for example Bregman iterations as proposed by Osher and coworkers [16]. To 
incorporate this iterative regularization method in our algorithm, we use their ” adding-back- 
the-noise” formulation such that the update for u in the previously introduced routine is 
replaced by 

u n+l = + h k + u n -TAy n+1 ). (30) 

Besides, the existing implementation is extended by an outer loop over a given number of 
Bregman iterations in which h is updated by h k+1 — h k + / — u after each complete cycle of 
the inner loop. 
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5 Computational Results 

In the following we present some results for the cases of image compression and denoising 
discussed above. As an example we chose the frequently used image ”Trui” (257 x 257 pixels), 
making the approach comparable to previous results such as m ■ However, since the size of 
this image does not correspond to modern HD resolutions, we also created two similar images 
with sizes of 1024 x 1024 and 4800 x 4800 pixels, respectively. 

5.1 Image Compression 




Figure 1: Comparison of the non-zero gradient ratio for different parameter values 

We start by discussing the compression of the cartoon part of an image, which is illustrated 
in Figure [l] for the Trui test image. We plot the relative error vs. the non-zero gradient ratio, 
which means the number of pixels with non-zero v divided by the total number of pixels. 
On the left we show a comparison of the sparse vector field (SVF) model with the classical 
ROF model, which demonstrates the improved compression properties. On the right we 
plot the results for the variational SVF model compared to the Bregman iteration (numbers 
correspond to Bregman iterations), which illustrates that no significant improvement can 
be obtained with respect to compression by the latter. In Figure [2] we display the results 
of the compression and the corresponding vector fields for A = 10. One observes that the 
support of v corresponds well to an edge indicator, confirming the relation to the approach in 
m • The reconstructed image seems to preserve the main edges well, but does not have the 
strict piecewise constant behaviour as total variation regularization, which seems attractive 
for further reconstruction tasks. 



Figure 2: From left to right: Original image, vector field v in x- and y-direction, norm of v, 
and the corresponding reconstruction for A = 10. 
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We also investigate the behaviour for higher resolution. In order to mimic increasing 
resolution we simply downscale the test images to r times the number of original pixels, 
r E (0,1]. We then perform compression at fixed error tolerance (corresponding to constant 
A when appropriately scaled) for the images of different size and finally plot the non-zero 
gradient ratio vs. r in Figure[3j Our expectation that due to continuum limit and the potential 
convergence towards a concentrated measure the ratio decreases with increasing resolution is 
well-confirmed for the Trui image as well as for a similar image at higher resolution. 



Figure 3: Comparison of the Non-Zero Gradient Ratios of the SVF (A = 10), the SVFBregman 
(A = 1 and 5 Bregman Iterations) and the ROF (A = 10) algorithms for the Trui test image 
(left) and an additional test image of size 1024x1024 pixels (right). The latter test image is 
displayed in the middle. 


Finally we display the result of the SVF model for image compression performed on a high 
definition image and compare it to a jpg image with the same compression rate in Figure [4j 
Indeed we achieve an improved PSNR with the SVF model. We also mention that several 
further compression steps on v can be carried out in an analogous way to [El, which will lead 
to highly improved rates, but is beyond the scope of this paper. 



Figure 4: From left to right: Original image with a resolution of 16.1068 bits per pixel, our 
reconstruction at a compression of 1.1892 bits per pixel (A = 10) with a PSNR value of 
34-0767 dB and the jpg image at the same compression rate with a PSNR value of 33.3214 
dB, corresponding difference images to original one in bottom row. 
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Figure 5: Noisy image (Gaussian noise, variance 0.05, left), ROF denoising result (A = 5, 
middle left), SVF denoising result (A = 4, middle right), norm of vector field v in SVF model 
(right). 


5.2 Denoising 

We finally give an outlook towards other tasks such as denoising with sparse vector fields. For 
this sake we compare the model with results of the classical ROF model and choose in both 
cases A such that the PSNR to the original image is maximized. We illustrate the result in 
Figure [5| which appears to be representative for all our tests. One observes that the reduced 
staircasing in the SVF model compared to the ROF model is less visible, which is due to 
point like artefacts that were not present without noise. This results in a lower PSNR than 
for the ROF model, which is consistent for all our tests. The reason for the artefacts is that 
v is too sparse in this case and does not encode the contours anymore. This can be seen in 
the last image of Figure [5j 


6 Conclusion 

We have introduced the SVF model for image compression motivated by diffusion inpainting 
and found several interesting connections to TV-type regularization methods. The SVF ap¬ 
proach leads to significantly sparser vector fields than the gradients of total variation, which 
appears quite useful for compression tasks, and seems not to suffer from staircasing artefacts, 
which appears attractive for other reconstruction tasks. However, the denoising performance 
of the SVF model is not convincing, since it creates point artefacts at reasonable choice of 
the regularization parameter. This is probably due to the fact that the norm induced by the 
corresponding regularization is too weak (there is an upper but no lower bound in terms of 
TV). For future improvement it seems natural to consider regularization on the curl of the 
vector field as well, such that v becomes again concentrated on contours rather than scattered 
points. In particular we suggest to study the more general problem 

^|| u — /Hi + || — w\\i+/3\\V x w\\i + 7||V • w\\i min (31) 

2 u,w 

for positive /? and 7 where again w = Vu — v (see Q). Due to the exact penalization 
properties of one-norms we expect that the ROF model corresponds to the case of /3 and 7 
sufficiently large, while the SVF model in this paper is /? = 0 and 7 sufficiently large. 
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